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Abstract 


FEASIBILITY OF A BLAST WAVE 
ATTENUATION STRUCTURE 


by Dale Richard Hartmann 


Chairman of the Supervisory Committee: Professor Ashley F. Emery 
Department of Mechanical Engineering 


This thesis begins with an overview of bombings in the United States, followed by the 
introduction of the Rankine-Hugoniot equations for blast wave pressure. The subsequent 
chapters develop the one-dimensional and two-dimensional Euler equations. These 
equations are the solved using the MacCormack finite difference algorithm. The basis of 
the investigation then begins by placing pole, shear plate and wedge obstacles in the path 
of the blast wave. The results of these simulations are interpreted and conclusions 
presented. Finally a synopsis of the existing results and cost analysis for structure 


hardening are presented. 
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PREFACE 


The inspiration for this thesis was developed over the course of two years. While at an 
assignment with the Defense Nuclear Agency, I was privileged to participate in several 
tests involving car bombs. I became very interested with the effects of the blast waves on 
different types of structures. After familiarizing myself with the current methods of 
protecting an existing structure from blast waves, | thought there must be a better way. 
This thesis is my investigation regarding whether there is a better way to protect existing 


structures. 
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INTRODUCTION 


BOMBS: THE UNKNOWN MENACE 


The threat of terrorist bombings is present every year. Although the United States has 
been relatively immune from terrorist bombings such as those in Ireland or Israel, 2577 
bombings occurred in the United States in 1995 alone. Bombings involving improvised 
explosive devices such as pipe bombs numbered 1,562 in 1995. The majority of the 
bombings in 1995 were small in size but produced $105 million damage and 937 
casualties (193 killed and 744 injured). The worst bombing in the history of the United 
States, the Murrah building in Oklahoma City, accounted for $100 million damage and 
786 casualties (168 killed and 518 injured). 


A cursory review of the above facts and events leads to the question: How can buildings 
that are possible terrorist bombing targets be protected from the effects of bomb blasts? 


That question 1s the basis of this thesis. 





CHAPTER 1: THE PROBLEM 


TYPES OF EXPLOSIVES 


Explosives classified as either: High explosives or Low explosives. High explosives 
possess certain characteristics that are vastly different than Low explosives. High 
explosives detonate, which means that when they are initiated a shock or blast wave will 
form. They also will burst or shatter materials near them, are capable of penetrating 


materials, and have the capability of lifting or moving objects. 


Low explosives, in contrast, do not detonate but burn very rapidly. Since Low explosives 
do not detonate, the pressure rise that is produced is usually smaller in amplitude but 
longer in duration than that of a High explosive. This combination tends not to produce 


an impulse type of blast wave but a slightly ‘softer’ shock to materials nearby. 


BLAST WAVE CHARACTERISTICS 


A blast wave generated as the result of the initiation of a contained High explosive is 
created in the following sequence of events. First hot gases with temperatures of the order 
3000 degrees C are generated. In concert with this temperature rise the pressure of these 
gases is of the order 100 to 300 kilobars. It is this second characteristic of high explosive 


initiation that is of concern here. 


This hot high-pressure gas then expands into the surrounding atmosphere. As expansion 
occurs, the air surrounding the expanding gas is forced outward. Since air is a 
compressible medium, a layer of air adjacent to the outer edge of the expanding gas is 
compressed. It is this layer of compressed air, which is called a blast wave. As the hot 
gas expands and cools, its pressure falls. The pressure of the layer of compressed air also 


falls as it is pushed farther outward from the point of detonation. As the gas continues to 





cool and expand, at some point in time and space the pressure will fall below atmospheric 
pressure due to the momentum of the layer of compressed air. This slight negative 
pressure is called overexpansion and results in a negative phase of the blast wave 
whereby the outward flow of gases and air is reversed. Eventually after some number of 


oscillations, equilibrium is reached and the motion of the air stops. 


The speed at which the blast wave travels is different depending on the medium in which 
it is traveling. For a bomb resting on the ground, two waves will exist: one traveling 
through the ground and another traveling through the air. The latter is what this thesis is 


concerned with. 


The magnitude of the impulse of an explosive blast wave is dependent upon many 


factors. Some of these factors are: 
Is the explosive cased? 
What type of explosive is used? 
The quantity of explosive used? 
What weather conditions are present at the time of detonation? 
What is the terrain? 


Are any structures nearby? 


PROTECTING EXISTING STRUCTURES 


THE TRADITIONAL APPROACH 
The traditional approach to protecting an existing structure from explosive blast damage 
is to harden it, that is to increase the structural strength of single or multiple components. 


Hardening of a building consists of several methods. 





The first and most rudimentary method is to harden the glazing, if any is present, on the 
exterior of the structure. Even a relatively small bomb can cause large amounts of injury 
and death without structurally damaging a building. This camage is achieved by 
shattering normal window glass and propelling the glass fragments at high speed 
throughout the interior of rooms along the outside perimeter of the building, anyone 
present in the rooms is impaled with large shards of glass. The use of tempered or 
tempered safety glass can minimize this effect. However, the use of blast curtains (heavy 
polymer curtains weighted along the lower edge) or the installation of Mylar film with 
reinforced window frames can eliminate this hazard. The problem is the curtains have to 
be kept closed at all times to provide protection, thus rendering the window useless as an 


architectural entity. In addition the Mylar is relatively expensive. 


The most complex, and best, method is to harden the entire structure against the effects of 
explosive blast waves. This generally consists of a major reconstruction of the interior 
structural elements of a building or encasing the building inside of another more robust 
exterior structure capable of withstanding the blast wave. This method, although highly 
effective, is also difficult to design, disruptive to the occupants, time consuming to 


construct and very expensive. 


THE NEW APPROACH 

The approach that this thesis investigates is the use of a blast wave attenuation structure. 
The form of the blast wave attenuation structure would be such that as the blast wave 
travels through it the energy of the wave is dissipated. Ideally, the energy of the wave 
would be lowered to such a degree as to produce and overpressure of no more than one 
half atmosphere. This lowered pressure, although still high enough to cause some 


damage, is safe for most commercial types of construction. 


The forms of the attenuation structure to be investigated include shear plates, a field of 
poles, and of wedges. The success of a type of structure will depend upon its ability to 


attenuate the pressure in less than 15 meters, constructability, low cost, and aesthetics. 





CHAPTER 2: THE QUESTION 


MODELING BLAST WAVES 


BLAST WAVE PARAMETERS 

The modeling of a blast wave must address two issues. The first is the modeling of the 
static aspects of the explosion. These static aspects include pressure, density, 
temperature and shock wave velocity. The second, and more difficult, is the modeling of 
the dynamics of the blast wave motion and interactions with objects. These dynamic 
aspects are the same quantities as the static aspects but are concerned with how these 


quantities change with time and position. 


STATIC ASPECTS 

It is important to model the static aspects of an explosion since this 1s what provides the 
driving force for the dynamic solutions and simulations. The first static aspect to be 
modeled is the pressure generated by the detonation of the explosive. The pressure 


generated is dependent upon several factors: 


Type of explosive 
Distance from explosive 
Position of explosive relative to the ground 


The type of explosive is important since each different type of explosive contains a 
different amount of energy per unit mass. These differences are summarized in the table 


below. 





Table 1 Explosive Equivalents 


Explosive Mass Specific Energy TNT Equivalent 
Qx(kJ/kg) (Qx/Qm7) 

RDX (Cyclonite) 5360 1.185 
HMxX 5680 1.256 
Nitroglycenn (Liquid) 6700 1.481 
TNT 4$20 1.000 
Blasting Gelatin 4520 1.000 
Nitroglycenn dynamite 2710 0.600 
Semtex 5660 20) 
Compound B 5190 1.148 


As can be seen from the table above, the method of normalizing the energies of different 
explosives relative to that contained in TNT has been developed and is universally 


accepted. 


The relationship between the range and the TNT equivalent charge mass is known as the 
scaled distance and 1s given by the following equation: 

=e (1) 
Where W, the universal symbol for TNT equivalents, denotes the mass of explosive and 
R denotes the range from the detonation of the explosive to the point of interest. Rankine 
and Hugoniot produced the first analytic solutions for blast wave front parameters in 
1870 for normal shocks in ideal gases. These solutions were later expanded by Brode to 
determine the peak static overpressure in the near and medium fields of the blast wave. 
The near field is the region where the overpressure is higher than 10 bar (10° Pa). The 
medium field is the region where the overpressure is less than 10 bar (10° Pa) and higher 


than 0.1 bar (10* Pa). The equation for the overpressure in the near field is given by 


6.7 
DiS & + 1} (100000) (Pa) (2) 


The equation for the overpressure in the medium field is given by 











- 455 9585 
ed eee 


eee - 0019) (100000) (Pa) (3) 


The absolute pressure of the detonation is given by 


Wee heetule a ee + aient (4) 


The equations developed by Rankine, Hugoniot and Brode, shown above, dealt with 
shock waves in free air. These shock waves were allowed to propagate in all directions, 
or, more simply, in a spherical manner. Since | am concerned with terrorist bombings, 
the majority of which involve explosives placed on or near the ground, a correction factor 
is needed ON Pabsotute to account for the reflection from the ground. If the ground were an 
ideal surface for reflection this correction factor would by 2. However, the ground will 
respond in two ways that are not ideal. The ground underneath the explosion will 
compress. The area around the explosion will undergo a fracturing process, which will 
result in particles being ejected into the air. There is no analytical method of quantifying 


these responses, but empirical evidence suggests a correction factor given by 


jb = SHE rite (S) 


This correction factor is applied to the pressure calculations for an explosion to provide 
the input pressure for all analytical and finite difference calculations. The Mach number 


of the detonation shock wave is given by 


| 1/2 
M, -|{(—2-- | fe ) (6) 
oer ay a I 


The units in the equations below are correct for pa expressed in Pa. 





The density of the gas behind the detonation shockwave is given by 





p> = ipo (Kg/m’) 
2 | 
Pern 
y+ M 
The velocity of the detonation shockwave is given by 


U., =< : [w, -+.]| (m/s) 
: Vaal M. 


Where c; 1s the ambient speed of sound 








1/2 
a = (; onlient (m/s) 


ee 


The temperature of the gas behind the detonation shockwave ts given by 





Pp 
Tf, =—~* (K 
=a (k) 
Where 
R 
ae 4 _(J/ke K 
uw. | g K) 
With 


R, = 8314.51 (J/kmole K) 


MW 


air 


= 29 (Kg/kmole) 


The speed of sound behind the detonation shockwave is given by 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 








1/2 
mL 
es -(«, . F = (m/s) (14) 


ambient 


DYNAMIC ASPECTS 


THE EULER EQUATIONS 


Air has a small thermal conductivity («) and also a small viscosity (v) therefore the terms 
in the Navier-Stokes equations that depend on «x and v can be neglected, which leads to 
equations in which the convective terms dominate and the fluid (in this case air) is treated 
as inviscid. The inviscid 2-D Navier-Stokes equations are called the Euler equations and 
are the equations dealt with in this paper. These equations are much simpler in form and 


also programming effort. 


The Euler equations are written as 


DF, _o (15) 
A Ee 
where 
p pu pv 
os puy 
pas pee Gan (16) 
py puy pv +p 
ie u(E + p) V(E + p) 


and 
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E= per hw 7) Gi) 


The Euler equations represent a system of equations that conserve mass, momentum and 
energy. These equations must be solved as a system to capture nonlinearities such as 


shocks. 


The Euler equations also need an equation of state, which 1s the ideal gas law. 


p= pk? =(y-1 B+ pw +v°)| (18) 


oe ™y 
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CHAPTER 3: ONE DIMENSIONAL EQUATIONS 


THE 1-D EULER EQUATIONS 


The first step in modeling the blast wave is to construct a one-dimensional model with 


time and distance as the parameters. The 1-D Euler equations are written as 


a Cx 
where 
p pi 
O=| pu f= pue + p (20) 
E u(E + p) 
and 
a e+’) (21) 
2 


The equation of state is 


p= pRT =(y- Ne - > Aw’) (22) 


This system of equations yields the following three equations 


Conservation of mass 


& , AP) _ 9 (23) 
ak 





Is 


Conservation of momentum 


Apu) i Apu a p) ae (24) 
ot OX 


Conservation of energy 


Eé. OU E + p) 
at Ox 


[| 
© 


(25) 


FINITE DIFFERENCING 

To evaluate this system of equations the method of finite differencing is used. Finite 
differencing 1s a powerful numerical technique that can be used to solve partial 
differential equations. Finite differencing algorithms come in two varieties: explicit and 
implicit. An explicit algorithm solves each equation in turn and then applies those results 
to the initial values for the next time step. An implicit algorithm solves the system of 
equations simultaneously via matrix manipulations. I chose the explicit form of finite 


differencing instead of the implicit form for the following reasons: 


Ease of coding 
Lower computational requirements 
Ability to run on desktop computing systems 


The first step in constructing a finite differencing program is to choose an algorithm. 
Many different algorithms exist but to be useful for the problem at hand the chosen 
algorithm must be able to “capture” the shock discontinuities well and also handle 


surface interactions. Several algorithms meet these requirements: 


Forward Time Centered Space (FTCS) 
MacCormack 


Harten-Ye 





Ist Order Roe 


| initially investigated the use of the FTCS algorithm, however I was not pleased with the 
shock “capturing” abilities of this algorithm. This algorithm smears or spreads the shock 
discontinuity over 4 to 5 Ax intervals, it also has dispersion errors which appear as 


oscillations at the shock discontinuity. The following figure demonstrates these points. 
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Figure 1 FTCS Output 


The next algonthm I investigated was MacCormack’s. This algorithm is a two-step 
method that is known for “capturing” shock discontinuities very well. In contrast to the 
FTCS algorithm, MacCormack’s does not have the dispersion errors at the shock 


discontinuity. The following figure demonstrates these points. 
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Figure 2 MacCormack Output 


The Harten-Yee and Ist Order Roe methods were abandoned early on due to their 
complexity and coding difficulties and because they have been shown not to be 


substantially superior to MacCormack’s method. 


Due to the simplicity and shock capturing abilities of the MacCormack algorithm, | have 


chosen to use it as the basis for the remainder of this thesis. 


The method of the MacCormack algorithm is contained in the two steps below: 


Predictor step: Q, = OF eae (26) 


Corrector step: on = 5(2; + Q,") — av F. (27) 


wed 


where 





IS 


Roa ah 


x jtl j 
Vee eel 
For the one dimensional case the MacCormack algorithm equations become 


Conservation of mass 


BD , Aet) _ 
a x 


— Nt : n 
Predictor: , =p; en (a es — (pu), 


x 








Conscion p= sens “| 


2 2Ax ey - pe.) 


Conservation of momentum 


Conservation of energy 


(28) 


(29) 


(30) 


(31) 


(32) 


(33) 


(34) 
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& OU E + p) =, (36) 


Predictor: iB, _ 5. = - 
Y 


At (pu) r+. , n n 
[Au"e,, iP js )-—+ 
P js) - 


n+] E. ae a AE ae oa a3 
OIG CLONE ie ee —(E; Sees Ds) (38) 
2 Die 
The equations for the conservation of momentum and energy above contain a pressure 
term that is evaluated using the equation of state. The equation of state is the ideal gas 


law and is given by 


p=(y- le er) (39) 


p 


In the above form the equation state can be evaluated using the results of the 


MacCormack algorithm. 


The final concern for any finite differencing algorithm is stability. The use of the 
Courant-Friedrichs-Lewy or CFL number is the most widely used method of controlling 
stability. The CFL number arises out of the results from a Von-Neumann stability 
analysis. The theory of Von-Neumann utilizes a Fourier transform to transform the finite 
difference solution into a wave space solution. The amplitudes of the waves in this 
resulting solution will grow or decay based upon the particular finite difference algorithm 
chosen and the value of the CFL number. The waves in the Von-Neumann solution that 
grow are unstable, those that decay are stable. The results of the Von-Neumann analysis 
require that the CFL number be less than or equal to one. The physical implication of the 
CFL number being less than or equal to one is that the sum of the amplitudes of the 


waves that decay is larger than the sum of the amplitude of the waves that grow. Thus if 





he 


the Courant number remains below the value of | the decaying waves will dominate and 


the algorithm will remain stable. 


AW 


aes (|x + c) At 


(40) 


where c is the speed of sound for the given pressure and density and CFL is the Courant- 
Friedrichs-Lewy number. To determine the time step used in the finite differencing code 
the CFL equation above is manipulated such that the following equation for a one- 
dimensional algorithm time step results. 


ene (41) 


(|u| + c) 


VALIDATING THE ALGORITHM 


At this point, after developing the equations for the one-dimensional case and 
constructing the finite difference algorithm, validation of the method is appropriate. To 
validate the code a model was constructed with a right traveling shock that reflects from 
an infinite wall at the right most-boundary. I selected this validation method due to its 


ease of coding and readily available analytic solution. 





Figure 3 Shocktube with right traveling wave 


The actual validation was done by allowing the finite difference solution to run long 


enough following the interaction with the right most wall and then comparing the 





L8 


pressure and density values from the model to the analytic values calculated from the 


analytic solution (Shapiro) 


The constant initial conditions in the program code are 


Table 2 Fixed initial conditions 


Temperature Gamma Density Pressure Speed of sound 


298 K 1.4 1.614 kg/m 101325 Pa 296 m/s 





The initial conditions input to the finite difference model were 


Table 3 Input initial conditions 


Courant Calculation Numberof Number of Distance Explosive Explosive 
number distance position time steps from charge mass 


steps 





The program calculated values behind the shock prior to the wall interface were 
determined using the Rankine-Hugoniot and Brode equations, equations 2, 3, 5, 6, 7, 8, 9, 


10 and 14. The values output from the program are tabulated below. 





Table 4 Calculated free stream shock values 


Static pressure of blast 5.7417x10° Pa 





Mach number of shock 2.2504 
Density 4.8729 kg/m 
Shock velocity 446.18 m/s 


Temperature 410.97 K 







Speed of sound 348.97 m/s 


The calculated values from the program behind the shock after the wall interface were 


Table 5 Calculated reflection shock values 


Static pressure of blast 2.246x10° Pa 


Density 12.179 kg/m 





To obtain the analytic values an iterative method was used. Following the reflection 
from the right most wall the shock is now traveling to the left. The analytic solution is to 


hold the shock stationary and iterate for the value of W. 
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Figure 4 Shocktube with stationary wave 


The equations required for the iterative solutions are 





wt u, 
Ci 42 
—_ (42) 
eee (7, ] (43) 
P2 ial 
ae (44) 
P2 1 2 i S 
Vay M 
2 
Ue a2 (wm, -4-| (45) 
vet | M, 


Using the equations above, equations 41 through 44, with the initial values listed in the 


table below, the results of the analytic solution will be those listed in the table below. 





a1 


Table 6 Analytic results for shocktube 


Initial values P . P2 (> 
446.18 m/s  348.1525m/s_ 574.17kPa 4.8729 kg/m*—s«d1.4 


Calculated values P3 (3 


180.5 m/s 2.0747 MPa 11.4963 kg/m” 





Comparing these analytic values with the finite difference calculated values yields an 
error of +6% for the density and +8% for the pressure. Which means that the finite 
difference program will yield a slightly conservative (higher than true values) result. 


This is acceptable for my purposes. 





CHAPTER 4: TWO-DIMENSIONAL EQUATIONS 


THE 2-D EULER EQUATIONS 


The addition of a second dimension to the 1-D Euler equations is fairly straightforward. 
The second dimension requires one additional equation to the matrix and two additional 


terms. The 2-D Euler equations are written 


Ze) ane & = 0 (46) 
a a Oh 
where 
p pu pV 
- +p PUuV 
o=|| ral"? | es] 7 (47) 
pv Puy OV a 
EB u( E a p) WE + p) 
and 


B= er t(u +’) (48) 


The equation of state is 


ek = (y — IE -= Alu’ 9) (49) 


This system of equations yields the following four equations 





Zo 


Conservation of mass 


—o a =0 (50) 


Conservation of x momentum 


Aon) dae +e) , Aa) 


5] 
Z (51) 
Conservation of y momentum 
Apu) Apw) Aer+r) (52) 
at Ox HY 
Conservation of energy 
At OWE + 
E AE+p) ME+D) _ te 


a Xx O 


The method of the MacCormack algorithm for two dimensions is similar to that for one 


dimension. The two dimensional algorithm is contained in the two steps below 


Predictor step: on =O, - Ai AF,” PO a, ") (54) 
Corrector step: On = ~(O,+ 2,.")-> a A (VF eG G;, | (55) 

where 
Nhe! (56) 





a ae oi (57) 





and 





> = Gia Gia 
Me = G; x = Ce 


For the two dimensional case the MacCormack algorithm equations become 


Conservation of mass 


Ae) Am) _ 
a ox OW 
a pu) , ~le,, (ev) le), 
Predictor: p,, = P;, — At (eu! Die Osa D9" 
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Corrector: pre = ae Pee = a PY ie nie 
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Conservation of x momentum 
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where 





p=(y- fe ie) . (66) 


and 
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Conservation of y momentum 
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where 
ye ee) 
p=(v fe : a (71) 
and 
p=(r- it Z leh (72) 
a 


Conservation of energy 
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The pressure terms in the conservation of energy equations above are evaluated by using 


the equation of state. The pressure terms are distinct for each direction. The x direction 


pressure term is given by 
p=(y+ | 2 {tad} (76) 


Similarly the y direction pressure term is given by 
2 
1{ (2) 
9=(v+4+1} E-— Td 
ae | p ) - 


The time step in two dimensions is a more complicated expression given by 


Nie = eC LL 
| Iv | l 
Woo eID sg ess 
Ax Ay ee 


(78) 
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where c is the speed of sound for the given density and pressure and CFL is the Courant- 


Friedrichs-Lewy number. 


VALIDATING THE ALGORITHM 


Applying reasoning similar to the one-dimensional case, an input pressure was applied to 
the x edge of the computational domain. The effect of this input pressure was the same 
as the calculations preformed by the one-dimensional code. Inputting the same 
parameters into both the 2-D and 1-D programs allowed comparison of the output 
matrices, which were identical. Therefore the 2-D code is validated since the 1-D code 
has been proved correct and acceptable. Following validation of the 2-D code for a linear 
pressure wave traveling in the y direction, the same procedure was applied for a wave 


traveling in the x direction. The same results were obtained. 


The results of this validation are shown in the tables below 


Table 7 1-D program calculated values 











Density 


4.8729 3.0026 Peal 1.6140 1.6140 


Momentum 


z782 743 Ee 


Energy 


ee 5co 1.0087e6 0.3144e6 0253960 0252526 






The results for a 2-D wave traveling in the x direction are given in the tables below 
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Table 8 2-D x-traveling wave program calculated 
values 


18729 roa 

18729 [era 
18729 [6140 
48729 610 
1879 6140 













Momentum 
x74 7a ae 
aT 7a a a 






EL a 

EL 2 

m2 [sms SOC 
Energy 











The results for a y traveling wave are in the tables below 
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Table 9 2-D y-traveling wave program calculated 
values 


1875 
50026 
7713 


1.6140 1.6140 1.6140 1.6140 1.6140 































 0067e6 
oa Tadet 
0253348 
0 2533e% 


974.5 974.5 974.5 974.5 


1.6140 1.6140 1.6140 1.6140 1.6140 
Momentum 
2742 2174.2 2174.2 2\G4 2 2174.2 
| 74 
753 
Energy 
ipo 566 | Bees ae 239566 19395e6 1.9395e6 


The results above were achieved with the following input conditions: 


Table 10 2-D validation initial conditions 


Time step x-grid Range to | Explosive Type | Explosive 
Charge Mass 


set 200 









CHAPTER 5: TESTING THE NEW APPROACH 


The validation of the 2-D code, presented in the previous chapter, was performed with 
linear wavefronts traveling in only one direction. While this situation might occur with a 
very large explosion at a great distance, in general this 1s an unrealistic case. Therefore, 
the use of symmetry was applied. The final step in the code development process was to 
place the explosion in one corner of the computational realm and allow the wave to 
expand in true 2-D fashion. Placement of the explosive in this corner allows for two 
planes of symmetry and reduces the number of calculations per time step by a factor of 
two. This reduction of calculations is also beneficial since it ignores the blast wave that 


travels away from the obstacles. 


The initial parameters for an explosion were placed into one corner of the computational 
realm and allowed to expand freely with no obstacles present. The results obtained were 
consistent with the expected outcome. This step was necessary to provide a test of the 2- 


D code in the true two dimensional environment. 


OBSTACLE BOUNDARY CONDITIONS 


To obtain meaningful results with obstacles in the path of the wave, a thorough 
understanding of the correct boundary conditions for the obstacles is needed. The 
minimum size of an obstacle is limited to a grid of 3 by 3 computational points. This 
minimum size is necessary to allow for the specification of conditions inside the obstacle. 
There is no limit on the maximum obstacle size, although it obviously cannot exceed the 


size of the computational realm. The figure below illustrates a minimum size obstacle. 
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Figure 5 Obstacle boundary conditions 


As can be seen in the figure above a 3 x 3 grid is the smallest obstacle possible to provide 
for the specification of quantities inside an obstacle. The figure above also shows what 
quantities exist at the points on the boundaries. The quantities that exist at the interior 


point are density and energy. 


For inviscid flow on rigid surfaces four distinct boundary conditions exist. The first is 
that the component of the x-momentum normal to the surface is zero. Similarly, the 
component of the y-momentum normal to the surface is also zero. The third condition is 


that the pressure gradient is zero. 
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The last boundary condition is that velocity component normal to the obstacle surface is 


ZEro. 


Another assumption needed to evaluate all parameters is constant energy inside the 


obstacle, #. The equation for energy is: 


= per a) (79) 


Since at this point the velocities are known, the energy and the pressure can be 
calculated. This equation along with the velocities allows for calculating p on the 


surface. 


The pressure at the surface of the body obtained by solving the ideal gas law equation, 


which is: 


Oe (y _ Ne ~= Alu’ +v')| (80) 
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RESULTS OF VARIOUS ATTENUATION STRUCTURES 


The following sections detail the results of running the simulation with poles, shear plates 


and wedges for obstacles. 


Before reviewing the results of the simulations it is helpful to examine the flow field 
results without any obstacles present. All of the results to be presented are at 30 time 
steps. The program does utilize variable time stepping so the results are not at the same 


total time. These results are shown in the figure below. 





Absolute Pressure 
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Figure 6 2-D flow field results without obstacles 


The figure above illustrates two points unique to numerical solutions. The first point is 
the gradient at the shockfront. The gradient of the shock is shown as the lines spaced 


closely together. These lines follow the curved path of the shockwave and travel at the 
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speed of the shock. The other point that is well illustrated, are the numerical oscillations 
inherent in this type of solution. The regions outlined behind the shockfront are, in fact, 


numerical oscillations of the solution and do not exist in the real explosion 


POLES 

The results for the field of poles were disappointing since it was the conceptual basis for 
this thesis. The pole obstacles, in the figure below, were set at 3 by 3 grid points with 
infinite height. The actual physical size of the obstacles varied with the total number of 
grid points and the calculation distance. For example a 10 by 10 grid, with a 4 meter 
calculation distance generates a pole that is 1.2 maters on a side. It is very important to 
ensure that the obstacle size is reasonable. In the above example to make the pole 
dimensions smaller a greater total number of grid points would be needed or a smaller 


calculation distance. 


Following the completion of many simulations the effect of poles on reducing the 
pressure of a blast wave it was found that the effect is minimal. The blast wave does 
compress the air on the blast side of the pole as I expected, but the shadow of the pole 
does not extend far enough in space to actually cause a decrease in the pressure of the 
wave after it has passed through the entire field. The cause of this phenomenon is the 
low viscosity of air. The calculations that are used for the simulation are inviscid, which 
allows the air to flow easily around a small object such as a pole. For the effect of an 
obstacle to be maintained farther downrange, the obstacle must be of a size large enough 
to produce a significant wake or a field of small poles set in a grid or staggered grid 
pattern that are close together. The field of poles would need to be set on the order of 
their diameter to be effective. Due to simulations running on a desktop computer system 
I could not generate a field of poles of sufficient density to demonstrate this directly. | 
am extrapolating a solution based upon the limited observations and conclusions | made 


following the simulations. The following figure illustrates this result. 
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Figure 7 Flow field results for pole obstacles 


This figure shows that the first pole, located at grid point (3,3) has a very large pressure 
gradient on its face. This gradient is produce by the compression of the air between the 
immovable pole and the moving air behind. This compressed air contains a great deal of 


energy that remains stagnant upon the face of this pole. 


The figure also shows that following passage of the shock through the pole field the 
gradient has been spread over a larger area, since the contour lines are farther apart, but 
the magnitude of the gradient has not changed. Therefore, the shock is not as steep as it 
was originally but is of essentially the same magnitude. Thus, the field of poles has 
‘softened’ the shock but has not attenuated the peak pressure. The interpretation of this 
result is that this spacing of poles, which is an example of many simulation runs, does not 


achieve the desired outcome. Extending this interpretation further leads to the conclusion 
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that no economical pole field would attenuate a shock to a level low enough to protect a 


Structure. 


SHEAR PLATES 

A shear plate is a thin long plate with the long axis ideally placed parallel to the direction 
of flow. Due to the uncertainty in the location of an explosive device, I placed the shear 
plates such that the length axis is perpendicular to the face of the protected building. The 
shear plates, in the figure below, are 3 grid points wide by 5 grid points long. | 
performed numerous simulations with other sizes and orientations but this combination of 


width, length and placement best demonstrates the effect of shear plates on the shock. 


The shear plates performed slightly better than the pole field largely due to their size and 
orientation. The first plate redirects the majority of the flow field in such a way that it no 
longer raises the pressure on the lower surfaces of subsequent obstacles. The pressure on 
the upper side of each shear plate is significantly less than the pressure on the lower side. 
This is the result of the lower plates acting as walls and shielding the upper plates from 
the blast pressure. The plates also redirect the flow by straightening it in the length 
direction of the plates. After the flow has passed the end of each plate it begins to 


diffract. The placement of the shear plates is critical to achieve this effect. 


The results of the simulation can be seen in the following figure. 
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Figure 8 Flow field results for shear plate : 
obstacles 


The figure above demonstrates the effect of shear plates on the flow. The lower plate 
redirects the majority of the flow along the lower edge of the computational realm and 
removes a great deal of the energy from the flow that impacts the middle and top plates. 
The flow then diffracts around the lower plate after it has passed. The gradient is again 
‘softened’, but its magnitude remains essentially the same as the original unimpeded 
shock. The large gradient in the lower left corner is produced by the shock impacting 
upon the face of the lower shear plate, which is immovable. This impact compresses the 
air and produces a large pressure peak in a very small area. The middle and top shear 


plates have similar gradients on their faces but the magnitude 1s much smaller. 


To provide significant protection for a building the location of the explosive would either 


need to be known or surmised before construction of the plates. The plates could then be 
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placed in such a way as to direct the flow away from the protected structure. 
Unfortunately, foresight is rarely this accurate and the construction of many plates (or 
walls for larger buildings) will be neither aesthetically pleasing nor particularly 


inexpensive. 


WEDGES 

A wedge is to two thin plates placed at an angle to each other and joined at the apex. The 
wedges used in this simulation consist of two plates 3 grid points by 5 grid points placed 
perpendicular to each other. The wedges were then placed into the computational realm 
such that the lower face of the wedge is diagonally across the path of the shock. The 


wedge in this position acts somewhat as a shield or wall in the path of the shock. 


The results for the wedge obstacles were similar to the shear plates. A large pressure 
gradient builds on the explosion side of the lower face and a shadow, or wake forms 
behind the wedge in the hollow interior. As the flow progresses past the end of the 
wedge diffraction begins and the shadow disappears. The main difference between the 
shear plates and the wedge obstacle is that, due to its size and geometry, the wedge does a 
better job of lowering the pressure in its wake. The wake, however, disappears rapidly, 
due to diffraction, after the shock passes the end of the obstacle. The orientation of the 
obstacle is not as critical as with the shear plates. The figure below illustrates these 


points 
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Figure 9 Flow field results for wedge obstacle 


The previous figure illustrates the effects of a wedge obstacle on the flow. The geometry 
of the wedge is such that one face is roughly perpendicular to the shock and a large 
pressure gradient builds on this face. This is shown by the closely spaced contour lines in 
the lower left corner between grid (4,4) and (7,2). The effects of this pressure gradient 
are that, the wedge absorbs a large amount of energy, the flow is split and a shadow or 
wake forms inside the faces of the wedge. The pressure gradient downstream of the 
wedge is ‘softened’, but its magnitude is, once again, similar to the unimpeded shock. 
Once again, following passage of the shock past the end of the wedge diffraction begins 


and the pressure begins to fill the shadow. This result is not unexpected. 


Since the orientation of the wedge is not as critical to its performance as that of the shear 
plates, wedges would be more a practical protection than the other obstacles investigated. 
However, due to the low viscosity of air the wedge would, in all likelihood have to be 


incorporated as an element on the exterior of the structure, and not a separate element. 
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This incorporation would lead to a structure with wedge shaped faces which, it is already 
known to demonstrate good resistance to surface loading from many sources including 
blasts, waves, running water etc. Finally this form of protection is neither inexpensive 


nor aesthetically pleasing. 
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CHAPTER 6: COST ANALYSIS 


The cost analyses for the obstacles investigated are moot. Since the pole obstacle does 
not work as hoped, unless many poles are placed closely together, in which case a wall 
would most likely be cheaper, | will present a synopsis of the available cost analyses for 
the conventional approach to blast protection. The following text is a compilation of 
others work and is not to be considered my own. These are my words but not my ideas or 


effort. 


Nuclear Disasters and the Built Environment contains an overview of the results of a 
testing completed in the 1950’s with above ground nuclear weapons. The results show 
that an overpressure of 5 psi (34.4 kPa) will completely demolish a conventional wood 
frame house and an overpressure of 1.7 psi (11.7 kPa) will result in serious damage. A 
later test series included strengthened wood frame houses. The first house, which was 
subjected to 4 psi (27.6 kPa) overpressure, sustained serious structural damage but 
remained standing. The second house was exposed to 2.6 psi (17.9 kPa) and suffered 
damage similar to the unreinforced house exposed to 1.7 psi. Also included in the second 
test were two brick houses. These were found to be of the same strength as the wood 
frame houses. One brick house was exposed to 5 psi (34.4 kPa) and sustained damage 
similar to the wood frame house. The second house received an overpressure of 1.7 psi 


(11.7 kPa) and sustained much less damage than the wood frame house. 


The strengthening of the houses for the test resulted in a 5% construction premium over 


unreinforced construction. 


The table below compiles the strengths of some common building materials. As can be 


seen the overpressure required for failure is quite low. 
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Table 11 Strength of common building materials 





Structural element Overpressure to cause failure 






(psi) 


Standard house construction 


Concrete wall panels 8”-12” thick 5-5-5 













Protecting Buildings from Bomb Damage contains an extensive and thorough cost 


analysis of hardening a new building to resist exterior explosions. I will summarize the 
results of the analysis. The model was based upon a 5% construction premium for 
hardening the building and also assumed a minimum 10% return on investment. The 
construction premium is the additional cost that would be incurred if the blast hardening 
were to be included in the construction of a building. The construction premium 
assumption is in keeping with the result from the tests completed in the 1950’s. The 
Output from the model was the lease premium based upon these assumptions. These 


results are shown in the table below. 
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Table 12 Construction-lease premium 


The lease premium in the table above is the additional amount that would have to be 







Assume construction premium 














charged lessees in a blast hardened building. This lease premium would recover the 
construction premium at the 10% return on investment. The table clearly shows that the 
resulting lease premium is less than the construction premium. Therefore it can be 
economically feasible to provide blast protection in new construction, if the developer or 


customer deems it necessary. 
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APPENDIX A: MACCORMACK 1-D CODE LISTING 


% Some constants needed for calculations (UNITS) 


%Courant number 
CFL=input('Enter the Courant number(O<CFL<1 0.5 is a good starting point): '); 
iene ||) 
error(‘The Courant number must be less than 1!') 
elseif CFL<=0. 
error(‘The Courant number must be larger than zero!') 


end 
Of ene ee 
%gamma for air 
g=1.4; 
eee, 
ioe area ee eee ed oe 
%rho = initial density 
(kg/m‘3) 
rho=1.614; 
er sn nae ae 
% Pin = initial pressure (Pa) 
Pin=101325; 
Dae rns Oe Se 
% T=initial temperature 
(K) 
T1=298; 
OE 
“calculate the speed of sound at the initial pressure 
(m/s) 
c=(1.4*Pin/rho)*(1/2); 
eee 
%h= distance over which to perform calculations (m) 
h=input (‘Enter the distance in meters over which to perform calculations: '); 
if h<=0. 
error(’You must enter a nonzero distance!') 
end 
a on een nee nne---------- 


“ostop= number of position steps per distance 
stop=input (‘Enter the total number of x grid points(101 is a good starting point): '); 
if stop<=0. 
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error(’You must enter a nonzero number of x grid points!') 
end; 


Yocalculate the x step size 
dx=h/(stop- 1); 


i 
% Print output of values on screen 
Uy SE eee 
fprintf('------------------------------------------------------------------ \n') 
fprintf(' Calculation values 
\n') 

Fennell Courant number 1S. ---....,....2200-c..0-......-.-....701.20\in' ,CFL) 
fprintf(‘The distance to calculate over 18........0..0c0ccc ce %3.2f m\n',h) 
fprintf(‘The number of x grid points iS...........0.0ccceee %4.0f\n',stop-1) 
DIMI Me X SUC ISIZE 1S. oe ea ete. stskeicecsa acess. sv % 1.4f m\n',dx) 
fprintf('------------------------------------------------------------------ \n') 
aera TT eg cee ae A ch AE A i Hee Ne Hole Ate ke Rk KE A EE 
*\n') 
fprintf(' Conditions ahead of the shock 
\n’) 
Maportobim Cli MGrATTAITIA 1S pee otter tees eadeecsovekee cesses .teess %1.1f\n',g) 
omni he initial pressure 1s,...........................%3.3¢ Pa\n',Pin) 
ipaiemen| Pie initial ENSity 18,,..........26..00080.2..c1-0.0-..0- % 1.4f kg/m‘%3\n' rho) 
fom e miielall temperate 1S 4... .cst.s.. de... %3.3g K\n',T1) 
foliar Whe speed Of SOUNd IS...))s ieee. %3.4f m/s\n',c) 
Hip Rewari aA AA AEE PR Re He ee Hee OE oko Seok A ok Ok 
*\n') 
(meio seis) 2d eo ne 
endt=input('Enter the number of time steps to calculate over: '); 

if endt<=0 

error(’You must enter a number of time steps larger than zero!') 

end 
a aa eo 


% define vector size 


oldp=[1:stop]*0.0; 
newp=| | :stop]*0.0; 
oldr=[1:stop]*0.0; 
newr=[1:stop]*0.0; 
barr=[1:stop]*0.0; 
oldru=[1:stop]*0.0; 
newru=| |:stop]*0.0; 
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barru=[1:stop]*0.0; 
newu=| I :stop]*0.0; 
oldE=[1:stop]*0.0; 
newE=[1:stop]*0.0; 
barE=[ 1:stop]*0.0; 
loc=[1:stop]*0.0; 


% Convert x grid points to actual distances 


for x=1:stop 


loc(x)=(x-1)*dx; 
end 
ee PI ee eee 
% Calculate Peak staic overpressure Ps based on charge size 
% and range to charge the driving force for the 
% calculations to follow 
i anu een se no 
% r=range to charge 

(m) 
r=input('Enter the range to the charge in meters: '); 

if r<=0 

error(’You must enter a range larger than zero!') 

end 
Ue 
% Print out table of explosive types 
% T=type of explosive 
fprintf('-------------------------------------------------------- \n') 
fprintf(’ Choose the explosive using the table below\t\n') 
fone OMAPOUNG: Bo. .o.c.ndet ees ves: Ln") 
IgoTi MMM ey CLONILE) 2220.02, o0eeee geese c.-.c.---2-.2\N') 
MifOOToNIAUOH  UReA VAD ING case e224 -2 neh ay sthavicsiseees sess setss. 3\n') 
Omni Natmogby Cent (MGUIG)....2-.-.00es-c......-<--.4\0') 
HifOMTalioat IM BORON Neste csr) OC .n msn ss Sasso 5\n') 
NOTA LMG es ASTUTE LAURIN ccc)... cates ca. 0cs0e00- 1. -2casce. 6\n') 
fprintf('60 percent Nitroglycerin dynamite ..............0...... 7\n') 
IiposIALIOM OR LIVE KOM sts cb4cc 2s.nndvadeacse-cvosersa certs 8\n') 
fprintf('-----------------------------—-------------------------- \n') 
0 ees aera 
% Ask for whiat type of explosive to use 
a ee Ds 
T=input(‘enter the type of explosive to use: '); 

if T==1 


S=1.148; 
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elseif T== 

S=1.185; 
elseif T==3 

S=1.256; 
elseif T== 

S=1.48 1, 
elseif T== 

S=1.000; 
elseif T== 

S=1.000; 
elseif T==7 

S=0.600; 
elseif T==8 

S=17250; 
else 
error(’You must enter an explosive type!’) 
end 


fprintf(‘The TNT conversion factor for this explosive is: %1.3f\n',S) 


% Ask what charge mass to use (kg of TNT) 
Uf 
M=input('Enter the mass of the charge in kilograms: '); 

if M<0O 

error(’You must enter a mass larger than zero!') 

end 
W=M*S; 
eee 8 _.___.__-------- 
% calculate peak static overpressure (Bar) 
z=r/(W“(1/3)); 
Ps=(6.7/(z*3)); 

if Ps<10 


Ps=((0.975/z)+(1.455/(2%2))+(5.85/(2%3)))-0.019; 
else Ps=(6.7/(z*3)); 


end 
if Ps<0.1] 
Ps=Ps; 
end 
a ee i ae 
% Convert Ps(the static overpressure) from Bar (Pa) 


% to Pascals from overpressure to absolute pressure 
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% and apply conversion factor for ground burst 


Sa a tee 
% Print output of values on screen 
fprintf('-----------------------------—------------------------------------ \n') 
fprintf(’ Explosive values 

\n' 
fprintf(‘The range to the charge from the x=0 point is..............%4.2f m\n',r) 
fprintf(‘The mass of the charge in TNT equivalents is.............%4.3f kg\n',W) 
fprintf('The static overpressure of the charge is................... %1.3g Pa\n',Ps-Pin) 
fprintf('------------------------------------------------------------------ \n') 
ie eee 
% Calculate intial velocity using Ps 
fy re wen ee 


Mach=sqrt(((Ps/Pin)-1)*((g+1)/(2*g))+1); 
Vin=c*(2/(g+1)*(Mach-(1/Mach))); 


% Print output of values on screen 
fprintf(‘The initial wave velocity I5..............0c:cce %4.4f m/s\n', Vin) 
fprintf(‘The shock wave Mach number Is...............0..0c0008 %2.4f\n',Mach) 
fprintf('------------------------------------------------------------------ \n') 

(tera eee eye 

% Calculate values behind the shock 
cae ee 


rho2=rho/(1-(2/(g+1))*(1-(1/Mach’2))); 
T2=Ps/(rho2*(8314.51/29)); 
e2m(G, 221 2) 1 1)00,5; 


Map Tol Me a ae at aM I ER A AR ok A A A A OE I Ok Hh ee EE 


*\n') 


fprintf(’ Conditions behind the shock 

\n') 

Hi OUUMNCL AE ME MES SUING MG hee eee ois. ss sence -vacnis ene: %1.3g Pa\n',Ps) 

iif CE(NS Gy 1S eee eee ec chev sxsswe %1.4f kg/m‘%3\n',rho2) 
Hemme ME TCTIRMERALNTC 1S yey... k eee icbsanacdhoonk. %3.3g K\n',T2) 
POUT MUNN SPEC OleSOUNG 158... cce6 es. 5..-.0-0000-0e0reov ee %3.4f m/s\n',c2) 


ae Mel Tecan Pactaatace fer Trem, hE Ae a a a Oe RE Hee ee oe Ae ee oe oe oA he A he 


*\n') 
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for x=1:stop 
oldp(x)=Pin; 
oldr(x)=rho; 
oldru(x)=0; 
oldE(x)=Pin/g 1; 
end 


% Define boundary conditions at grid point x=1 

ee a ne we ee === 
oldp(1)=Ps; 
oldr(1)=rho/(1-(2/(g+1))*(1-(1/Mach*2))); 
oldru(1)=oldr(1)* Vin; 
oldE(1)=(Ps/g1)+(0.5*oldru(1)*2/oldr(1)); 


1) 

% Calculate the new speed of sound based on new pressure 

% and density values 

SINE asin Pee OD 
newc=sqrt(1.4*oldp(x)/oldr(x)); 

fs La eo 

% Calculate the initial time step 

1 lel ley ee 

dt=CFL*dx/(newc+toldru(1)/oldr(1)); 

fom che imital time StepiS.......2--.....--.....-- %1.5e s\n',dt) 

Di ee ee eee ee 

% Initialize total time counter 

5 genre es eon 

tt=0; 

fe ee 

% Define predictor values at x equal 1 

poe eS ee ee 
barr(1)=oldr(1 )-(dt)*(oldru(2)-oldru(1)); 
barru(1)=oldru(1 )-(dt)*(((oldru(2)*2/oldr(2))-... 

(oldru(1 )*2/oldr(1)))+(oldp(2)-oldp(1))); 
barE(1)=oldE(1)-(dt)*((oldE(2)+oldp(2))*oldru(2)/... 
oldr(2)-(oldE(1)+oldp(1))*oldru(1)/oldr(1) ); 

0 un ee os nec e 

% Adjust x grid end point for finite differencing 

stop=stop- 1; 
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% Finite difference for density, momentum and energy 


for t=1:endt 
dtmin=1.; 
CFLmin=CFL; 
dtdx=dt/(dx); 
%Calculate total time 


tt=tt+dt; 
for x=2:stop 
apy Sepa ee ee 
% PREDICTOR 
fi i a 
% Solve as a function of time and position for increasing 
% time 
RN Sree Salt Snare Oe Se tS 
% Density(r) 
barr(x)=oldr(x)-(dtdx)*(oldru(x+1)-oldru(x)); 
Ye op se tg 
% Momentum(ru) 
barru(x)=oldru(x)-(dtdx)*((oldru(x+1)%2/oldr(x+1))-... 
(oldru(x)*2/oldr(x))+(oldp(x+ 1)-oldp(x))); 
Oe A 
% Energy(E) 
barE(x)=oldE(x)-(dtdx)*((oldE(x+1)+oldp(x+1))*... 
(oldru(x+1)/oldr(x+1))-(oldE(x)+oldp(x))*... 
(oldru(x)/oldr(x))); 
Sennen rey Sage 
% CORRECTOR 
anon ee NE i 
% Density(r) 
newr(x)=(barr(x)+oldr(x))/2-(dtdx/2)*(barru(x)-barru(x-1)); 
% Left wall boundary conditions 


newr(stop+ 1 )=newr(stop-1); 
newr(stop)=newr(stop-1); 
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En 2) te 
% Momentum(ru) 
newru(x)=(barru(x)+oldru(x))/2-(dtdx/2)*((barru(x)*2/... 
barr(x))-(barru(x-1)*2/barr(x-1))+(g 1 *(barE(x)-... 
(.5*(barru(x)*2/barr(x)))))-(g 1 *(barE(x-1)-(.5*... 
(barru(x-1)*2/barr(x-1)))))); 
% Left wall boundary conditions 
newru(stop)=0; 
newru(stop+1)=0; 
Ey mercenaria 
% Energy(e) 
newE(x)=(barE(x)+oldE(x))/2-(dtdx/2)*((barru(x)/... 
barr(x))*(barE(x)+(g 1 *(barE(x)-(.5*(barru(x)%27/... 
barr(x))))))-((barru(x- 1)/barr(x-1))*(barE(x-1)+... 
(g 1 *(barE(x-1)-(.5*(barru(x-1)*2/barr(x-1)))))))); 
Se Rrra eee an 
% Calculate pressure using the continuity equation for the 
% ideal gas law and the new values for energy, momentum 
% and density 
newp(x)=g1 *(newE(x)-(.5*(newru(x)*2/newr(x)))); 
newp(stopt+1)=newp(stop-1); 
Wi ie 
% Calculate velocity using the new values for momentum 
% and density 
newu(x)=newru(x)/newr(x); 
pom alae 
% Check stability 
% Calculate the new speed of sound 


newc=sqit(g*newp(x)/newr(x)); 





% 


% ee ee ens 


% 


% =. - ee 


% 
% 


dt=dtmin; 
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Calculate the new time step dt 


dtnew=dx/(newu(x)t+newc); 
if dtnew < dtmin 

dtmin=dtnew; 
end 
CFLnew=(dtmin/dx)*((newctoldru(x)/oldr(x))); 
if CFLnew<CFLmin 

CFLmin=CFLnew; 


Return for next x 
end 


ee ome ee ee ee ee a SS SSS SS SS SSS SSS SS SS SCS SSeS] =——-— 


Apply numerical boundary scheme to calculate all values 
at x=1 and x=stopt1 


newr(1)=oldr(1); 

newE(1)=oldE(1); 

newru(1)=oldru(1); 

newu(1)=oldru(1)/oldr(1); 

newp(1)=oldp(1); 

barr(1)=oldr(1)-(dt)*(oldru(2)-oldru(1)); 

barru(1)=oldru(1)-(dt)*((oldru(2)*2/oldr(2)-... 
oldru(1)*2/oldr(1))+(oldp(2)-oldp(2))); 

barE(1)=oldE(1)-(dt)*((oldE(2)+oldp(2))*oldru(2)/... 
oldr(2) -(oldE(1)+oldp(1))*oldru(1)/oldr(1) ); 

newE(stop+1)=oldE(stopt1); 

newu(stop+1)=newru(stop+1)/newr(stopt1); 


Update time step and CFL(Courant number) 


CFL=CFLmin; 


fprintf(('CFL= %1 .3g\tdelta t= %1.5e s\tTotal time= %1.5e s\n',CFL,dt, tt) 





a) 


% Output the values for pressure, density and velocity 
% using graphs 

un enenereenareeeets Pies 55 kn 

% Pressure 


subplot(3,1,1), plot(loc,newp), ylabel(‘Absolute Pressure (Pa)’) 
% a ae 

% Density 

subplot(3,1,2), plot(loc,newr), ylabel('Density (kg/m‘%3)') 

% en ewww 

% Velocity 


subplot(3,1,3), plot(loc,newu), ylabel('Velocity (m/s)').... 
xlabel('Distance (m)') 


%Replace old values with new values and begin next time step 


oldp=newp; 
oldr=newr; 
oldru=newru; 
oldE=newE; 





APPENDIX B: MACCORMACK QUASI] 2-D CODE LISTING 
clear all 


% Some constants needed for calculations (UNITS) 


% Courant number 
CFL=input(‘'Enter the Courant number(0O<CFL<1 0.5 is a good starting point): '); 
beh ie lk 
error('The Courant number must be less than 11’) 
elseif CFL<=0. 
error(‘The courant number must be larger than zero!) 


end 
OMEN ee eee rise 0 oe oo. cols ons. ---u- 
“gamma for air 
p=1.4; 
elSeas 
Uy NE eee ie Se eer ak 
%rho = initial density 
(kg/m‘3) 
rho=1.614; 
0 Nn renee Pe r= nee re Sa 
%Pin = initial pressure (Pa) 
Pin=101325; 
Of ee 
%T=initial temperature 
(K) 
T1=298; 
0 
“calculate the speed of sound at the initial pressure (m/s) 


c=(1.4*Pin/rho)*(1/2); 


a Ie neces ene 
%h= distance over which to perform calculations (m) 
h=input (‘Enter the distance in meters over which to perform calculations: '); 


if h<=0. 
error('You must enter a nonzero distance!’) 
end 
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%stop= number of position steps per distance 
stop=input (‘Enter the total number of grid points in the x,y directions(31 is a good 
starting point): '); 

if stop<=0. 

error(’You must enter a nonzero number of x grid points!’) 

end; 


% calculate the x and y step size 
dx=h/(stop-1); 
dy=h/(stop-1); 


so ene Curse re ee 
% Print output of values on screen 
fprintf('------------------------------------------------------ \n') 
fprintf(’ Calculation values 

\n') 
fprintf(‘The Courant number is.. eee. Vol Jin Cr) 
fprintf(‘The x distance to calculate c ON Gti Seimei tia %3.2f m\n',h) 
fprintf(‘The y distance to calculate over IS...............c.cceee %3.2f m\n’,h) 
iPnmleheaiie NUMber Of X STIG POINTS 1S. ......-...2..0c2.00.22-2--- %4.0f\n',stop-1) 
HoMmimainemumbenr Of y Std POMS 1S. ......-.....0......<1--20-- %4.0f\n',stop- 1) 
WPUMIMMMENNG XESLEP SIZE 1S) ....4.6.c.av-c--,c0ueseteretbeers sss. %1.4f m\n',dx) 
IoMANMUMMelG Ve SLEP SIZE 1S. .c2.0c7.clees.d-f-seeec screens sees % 1.4f m\n',dx) 
fprintf('------------------------------------------------------------------ ') 
IVT pen UU acca gee fea Canaan ee ocr crane NC ME AR AE AC A eae a oe Hc co Hk oe a eee 
*\n') 
fprintf(’ Conditions ahead of the shock 
\n') 
Muestal Maia One MIAN Tcl Setanta omen nr ear lakiktechls vcs %1.1f\n',g) 
portant pall nie aime Tal NESSUNe 15.202. eres. -..cc...-..s. 2 %3.3¢ Pa\n',Pin) 
Mnf MeMMUILN US INENG sMMMMENAll (GONISILY 15.0.52.)200efeiase.-0seed0ese0 0s +s. %1.4f kg/m%3\n',rho) 
MO MMeeNemettal TeMPeratune 1S... 8 cn... .esie.t--...- %3.3g K\n',T1) 
MUMIMUM Me SPEC ON GOUMG IS.2..20000...ic.cc2-s..c02.5-200- %3.4f m/s\n',c) 


fiesta AAT cea aac Msi dic acre ci A Me rk ooh aI SA ok ah oh ke EE EA 


*\n') 
% ee oem nae enon aaannmaaeeeoaeos 


endt=input('Enter the number of time steps to calculate over: '), 
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if endt<=0 
error('You must enter a number of time steps larger than zero!’) 
end 

i aa eR Se Eo weeded 

% define vector size 

for x=1:stop 

for y=I:stop 


op(x,y)=0.0; 
np(x,y)=0.0; 
Olr(x,y)=0.0; 
Bicy 0.0; 
br(x,y)=0.0; 
Olru(x,y)=0.0; 
nru(x,y)=0.0; 
Olrv(x,y)=0.0; 
nrv(x,y)=0.0; 
bru(x,y)=0.0; 
mul x,y )=0.0; 
brv(x,y)=0.0; 
nv(x,y)=0.0; 
BE( x.y )—0.0: 
nE(x,y)=0.0; 
bE(x,y)=0.0; 


end 
end 


% Calculate Peak static overpressure Ps based on charge 
% size and range to charge the driving force for the 

% calculations to follow 

eee iarge otewuretes Siacs8 20 o ovetnoe sk 

% r=range to charge 


r=input('Enter the range to the charge in meters: '); 
if r<=0 
error('You must enter a range larger than zero!') 
end 
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% Print out table of explosive types 
% T=type of explosive 
fprintf('-------------------------------------------------------- \n') 
fprintf(’ Choose the explosive using the table below\t\n') 
Hips THIMTILAG COL IPS OUIINGIM Ey 40 go ee 222.0522 4.-¢4e0s ae eseeseoseeee 1 MM’) 
Ig OAR IMTEN (YO aCe VCNOMILE) 0.0. .occ ede -ss4.c2e.-+s-08000.2\N') 
Asfestol PH WEE LVL Neer Shae’. ne2 vines; cksanaxdescedatecensssiacles: 3\n') 
iotanmant INGtrOe ly Cem (MGUIC),...........cceccc<..ee seca 4\n') 
Ifo TaMMGN MUN lh teen... 2-2 se desea esteNlis iiss. 5\n') 
Rerun (een ASePAO Greatly 2... eaii.ilitie.d sits... cteac.. 6\n') 
fprintf('60 percent Nitroglycerin dynamite .................... 7\n') 
IIPS MINUS Xetra. cc cge ree ae STs xe clsens sx. 8\n') 
fprintf('-------------------------------------------------------- \n’) 
ae se 
% Ask for what type of explosive to use 
NTE eee ce 
T=input(‘enter the type of explosive to use: '); 
if T== 
S=1.148; 
elseif T==2 
S=1.185; 
elseif T==3 
S=1.256; 
elseif T== 
S=1.481; 
elseif T==5 
S=1,000; 
elsaif T== 
S=1.000; 
elo ail 
S=0.600; 
elseif T==8 
S=1.250; 
else 
error('You must enter an explosive type!') 
end 


fprintf(‘The TNT conversion factor for this explosive is: %1.3f\n',S) 
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% Ask for mass of charge 
(kg) 
ee wae a og esa ae 
M=input('Enter the mass of the charge in kilograms: '); 
if M<O 
error(’You must enter a mass larger than zero!') 
end 
W=Mi%S; 
ee 
% Calculate peak static overpressure (B) 
6 um en SPs 
z=1/(W”(1/3)); 
Ps=(6.7/(z%3)); 
if Ps<10 
Ps=((0.975/z)+(1.455/(z*2))+(5.85/(z%3)))-0.019; 
else Ps=(6.7/(z*3)); 
end 
goes 081 
Ps=Ps; 
end 
0 EE 
% Convert Ps(the static overpressure) from Bar (Pa) 
% to Pascals from overpressure to absolute 
% pressure and apply conversion factor for 
% ground burst 
GY ue eu lenses ee 
Ps=1.8*((Ps*100000)+Pin); 
Of 
% Print output of values on screen 
ee ae eg 
fprintf('-------——----—------------__----_--_--__--__--___-------.------ \n') 
fprintf(’ Explosive values 
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fprintf(‘The mass of the charge in TNT equivalents is.............%4.3f kg\n', W) 
fprintf(‘The static overpressure of the charge is...........0.00.... %l. 3g Pa\n',Ps-Pin) 
fprintf('------------------------------------------------------------------ \n') 

a oes 

% Calculate intial velocity using Ps 
0 cig ge ee Se 


Mach=sart(((Ps/Pin)-1)*((g+1)/(2*g))+1); 
Vin=c*(2/(g+1)*(Mach-(1/Mach))); 
% Print output of values on screen 


fprintf('‘The initial wave velocity is...........0.0.ce ee %4.4f m/s\n', Vin) 
fprintf(‘The shock wave Mach number is....................:00.....%2.4f\n',Mach) 
fprintf('------------------------------------------------------------------ \n') 

on 22 ee Ret nn 

% Calculate values behind the shock 
ef Nee 


rho2=rho/(1-(2/(g+1))*(1-(1/Mach’2))); 
T2=Ps/(rho2*(83 14.51/29); 
c2=(c*2*T2/T1)°0.5; 


MppestaV NH Te a as a al ae PE i neh ot oh FG oc i eae SE Ee He ee He ak a 2 2 Co 2 2k Ak 


=n) 
fprintf(’ Conditions behind the shock 
\n') 
PTI s RE SSVe aS Wye teneroeescecs see %1.3g Pa\n',Ps) 
Mpa mITMUA (GMA CMSIUY 1S hc. s ces satag sees eececseccoveses acer. %1.4fkg/m%3\n',rho2) 
IGIM MUU IDeTEMAMCTALULE 15,......,...00ee0rss0c2c0+000:-400eeeeeteceeee %3.3g K\n',T2) 
fprintf(‘The speed of sound is.. ee: ..%3.4f m/s\n',c2) 
Pn GLEEEEEe EEE EEE EEEEEELEE EEE EEE TEER LEER EE CELE LEER EEL 
*\n') 
Cen ee ee 
i / ei Seo este 
% Define initial conditions 
De ey 
for x=1:stop 
for y=I:stop 
op(x,y)=Pin; 
Olr(x,y)=rho; 


Olru(x,y)=0; 
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Olrv(x,y)=0; 
oE(x,y)=Pin/g1; 
end 
end 
% Ew ee Ue eee 
% Define boundary conditions at grid point x=1 y=1 
% ar a a 
y=; 
for x=1:stop 
op(x,y)=Ps, 
Olr(x,y)=rho/(1-(2/(g+1))*(1-(1/Mach’2))); 
Olru(x,y)=0; 
Olrv(x,y)=Olr(x,y)* Vin; 
oE(x,y)=(Ps/g1)+(0.5*(Olru(x,y)*2+Olrv(x,y)*2)/... 
Olr(x,y)); 
end 
% ah Pe Re eA R a a  B 
% Calculate the new speed of sound based on new pressure 
% and density values 
% a ng nea a ee cams ee Se 
ne=sqrt(g*op(x,y)/Olr(x,y)); 
% Nee eee ee eee eee eae eee 
% Calculate the initial time step 
% ee eee eee ee ee eee oa ene conan awoeand 


dt=CFL*(1/(abs(Olru(1,1)/Olr(1,1))/dx+abs(Olrv(1,1)/Olr(1,1))/... 
dy+ne*(1/dx*2+1/dy*2)*(0.5))); 


MO MIPENe ie wMItAlMNINLe SLED 1S f1..4ccct cscs. esacean eens %1.5e s\n',dt) 

1) ae NEINENIIRIRE cscs 

% Initialize total time counter 
uN OURS ene 
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stop=stop- 1; 


for t=1:endt 


% ee a a ee ee eee ee 
% NBS for column | and row 1 
Vip ee Soe eee 
dtmin=1.; 
dtdx=dt/dx; 
dtdy=dt/dy, 
ee a eae are Pepe eee eames wee 
% Calculate total time 
Od en ee Oe re ee 
tt=tt+dt; 
Se ee ee eee eS. 3. +-~-~= 
Non = |esuteye 
for y=1:stop 
eee en et eee eae cen- se ~-~ 
% PREDICTOR 
eae ee ee oo cca a none 
% Solve as a function of time and position for 


% increasing time 





pa 
% Density(r) 
Uf ee 
br(x,y)=Olr(x,y)-... 
(dt*(((Olru(x+1,y)-Olru(x,y))/dx)+... 
((Olrv(x,y+1)-Olrv(x,y))/dy))); 
br(x,stopt1)=br(x,stop); 
br(stopt 1,y)=br(stop,y); 
br(stop+1 ,stop+ 1 )=br(stop,stop); 
Of, SE es ie RON co arises Ae 
% Momentum(ru) 
serene ek pitt eno Meee it 
bru(x,y)=Olru(x,y)-dt*... 
((1/dx*... 
(Olru(x+1,y)*2/Olr(x+ 1,y)-... 
Olru(x,y)*2 /Olr(x,y) +... 
(op(xt+1,y)-op(xy) __)))t... 
lidy 
(Olrv(x,y+1)*Olru(x,y+1)/Olr(x,y+1)-... 
Olrv(x,y) *Olru(x,y) /Olr(x,y) _))); 
bru(x,stopt1)=bru(x,stop); 
bru(stop+1,y)=bru(stop,y); 
bru(stop+1,stop+1)=bru(stop,stop); 
A err enn nue 
% Momentum(rv) 
0) we 


brv(x,y)=Olrv(x,y)-dt*... 
(Citabee 
(Olru(x+1,y)*Olrv(x+1 ,y)/Olr(x+1,y)-... 
Olru(x,y) *Olrv(x,y) /Olr(x,y)_ ))+... 
(idye 
(Olrv(x,y+1)*2/Olr(x,y+1)-... 
Olrv(x,y)*2 /Olr(x.y)--... 
(op(x.y+1)-op(x,y) —_)))); 


brv(stop+1,y)=brv(stop,y); 
brv(x,stop+1 )=brv(x,stop); 
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brv(stop+1 ,stop+1)=brv(stop,stop); 


% a ee eee ee moon 


of, Energy(e) 


bE(x,y)=o0E(x,y)-dt*... 
((1/dx*... 
((oE(x+1,y)+op(x+1,y))*Olru(x+1,y)/Olr(x+1,y) -... 
(oE(x,y)top(x,y) )*Olru(x,y) /Olr(x,y)_)+... 
dy... 
((oE(x,y+1)+op(x,y+1))*Olrv(x,yt+1)/Olr(x,yt1) =... 
(OE(x,y)top(x,y) )*Olrv(x,y) /Olr(x,y)  ))); 


bE(stop+1,y)=bE(stop.,y); 


bE(x,stop+1)=bE(x,stop); 
bE(stop+1,stop+1)=bE(stop,stop); 


% return for next x 


end 
% return for next y 
end 
0% ee eG en a ee ee i a eee ee ee eee ee 
nr(:,1)=Olr(,1); 
nru(:,1)=Olru(:, 1); 
nrv(:,1)=Olrv(:,1); 
nE(:,1)=oE(,1); 
for x=2:stop+l 
for y=2:stop+l 
% ee ee a 
% CORRECTOR 
% 2 ae ee en ewe eee 
A Density(r) 
0% ee a a ee YS 


nr(x,y)=(br(x,y)+Olr(x,y))/2-... 
( dtdx/2*(bru(x,y)-bru(x-1,y))+... 
dtdy/2*(brv(x,y)-brv(x,y-1))); 
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nr(1,y)=nr(2,y); 


% Left wall boundary conditions 


% Momentum(ru) 


Px=g1*(bE(x,y)-0.5*((bru(x, y)*2/br(x,y))+(brv(x,y)*2/br(x,y)))); 
Pxm=g 1 *(bE(x- l,y)-0.5*((bru(x-1,y)*2/br(x-1,y))+(brv(x- 1 ,y)*2/br(x- 1,y)))); 


nru(x,y)=(bru(x,y)+Olru(x,y))/2-... 
((@idxi2 >=. 
(( bru(x,y) *bru(x,y) /br(x,y) +Px )-... 
( bru(x-1,y)*bru(x-1,y)/br(x-1,y)+Pxm)))+... 
(dtdy/2*... 
( bru(x,y) *brv(x,y) /br(x,y)-... 
bru(x,y-1)*brv(x,y-1)/br(x,y-1) _))); 


nru(1,y)=nru(2,y); 


% Momentum(rv) 


Daa ene eee Be ee Me Pee ae ee 


Py=g1*(bE(x,y)-0.5*((bru(x,y)*2/br(x,y))+(brv(x,y)*2/br(x,y)))); 
Pym=g 1 *(bE(x,y-1)-0.5*((bru(x,y-1)*2/br(x,y-1))+(brv(x,y-1)*2/br(x,y-1)))); 


nrv(x,y)=(brv(x,y)+Olrv(x,y))/2-... 
((decdag 27. 
( bru(x,y) *brv(x,y) /br(x,y) — -... 
bru(x-1,y)*brv(x-1,y)/br(x-1,y) ))+... 
(etd 25, 
( (brv(x,y)*brv(x,y)/br(x,y)+ = Py )-... 
(brv(x,y-1)*brv(x,y-1)/br(x,y-1)+Pym)))); 


nrv(1,y)=nrv(2,y); 
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% Left wall boundary conditions 


% Energy(E) 


nE(x,y)=(bE(x,y)+oE(x,y))/2-... 
((a@tdx 2". 
((bru(x,y) /br(x,y) *(bE(x,y) +Px )) -... 
(bru(x-1,y)/br(x-1,y)*(bE(x-1,y)+Pxm)))) +... 
(dtdy/2*... 
((brv(x,y) /br(x,y) *(bE(x,y)+Py )) -... 
(brv(x,y-1)/br(x,y-1)*(bEQ,y-1)+Pym)))) _); 


ee nares e ete J. 
0 ot ee 
% Right wall boundary conditions 
1p ne 
nr(x,stop+1)=nr(x,stop-1); 
% nr(x,stop)=n1r(x,stop-1); 
nrv(x,stop)=0; 
nrv(x,stopt+1)=0; 
nE(x,stopt+1)=oE(x,stop+1); 
nE(1,y)=nE(2,y); 
Uo 
“return for next x 
end 
Ree i 
“return for next y 
end 
te 
for Va estopa | 
for Kea leStOpan! 
0) a ces eee 
% Calculate pressure using the continuity equation for 
% the ideal gas law and the new values for energy, 
% momentum and density 
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np(x,y)=(g1)*(nE(x,y)-(.5 *((nru(x,y)*2/... 
nr(x,y))+(nrv(x,y)*2/n1(x,y))))); 


np(x,stop+1)=np(x,stop-1); 


0 re ene ree ee 2 
% Calculate velocity using the new values for momentum 
% and density 

Sees sn eo ee 


nu(x,y)=nru(x,y)/nr(x,y); 
nv(x,y)=nrv(x,y)/nr(x,y); 


ee a ee eS hee eee 


% Check stability 
% Calculate the new speed of sound 


dtn=CFL*(1/(abs(nu(x,y))/dxt+abs(nv(x,y))/dyt... 
ne*(1/dx*2+ 1/dy*2)*(0.5))); 
if dtn < dtmin 


dtmin=dtn; 
% end dtmin if statement 
end 
ea nee ena Mnenn ee eo et 
Oe 


YE ce pg ay 


%return for next x 
end 


“return for next y 
end 
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dt=dtmin; 


fprintf(‘delta t= %1.5e s\tTotal time= %1.5e s\n',dt,tt) 


% Output the values for pressure, density and velocity 
% using graphs 


% Pressure 


colormap(cool) 

surf(np) 

hold on 
pcolor((np)-Pin) 
title(’Absolute Pressure’) 
xlabel('X grid points’) 
ylabel(‘Y grid points’) 
Zlabel('Pascals') 


pause(.01) 
hold off 


% Replace old values with new values and begin 
% next time step 


Op—Np, 
Olr=nr; 
Olru=nru; 
Olrv=nrv; 
oE=nE; 





v2 


% Return for next t 
end 





APPENDIX C: MACCORMACK 2-D CODE LISTING 


clear all 


% Some constants needed for calculations (UNITS) 


% Courant number 
CFL=input('Enter the Courant number(O<CFL<1 0.5 is a good starting point): '); 
Hee ne 
error(‘The Courant number must be less than 1!’) 
elsetf CFL<—0. 
error('The courant number must be larger than zero!') 


end 
ae cee a nnn 
“gamma for air 
g=1.4; 
geo: 
os ho ea 
%rho = initial density 
(kg/m*3) 
rho=1.614; 
fe 
% Pin = initial pressure (Pa) 
Pin=101325; 
(mmm sen 2 
% T=initial temperature 
(K) 
T1=298; 
Up 
“calculate the speed of sound at the initial pressure (m/s) 


c=(1.4*Pin/rho)*(1/2); 


pn 
%h= distance over which to perform calculations (m) 
h=input (‘Enter the distance in meters over which to perform calculations: '); 


if h<=0. 
error('You must enter a nonzero distance!') 
end 
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Y%stop= number of position steps per distance 
stop=input (‘Enter the total number of grid points in the x,y directions(31 is a good 
starting point): '); 

if stop<=0. 

error('You must enter a nonzero number of x grid points!') 

end; 


“calculate the x and y step size 
dx=h/(stop-1); 
dy=h/(stop- 1); 


benno er ee 8 
% Print output of values on screen 
fprintf('------------------------------------------------------ \n') 
fprintf(’ Calculation values 

\n' 
in Mtinetd MMMe M@OUTANt NUMDEL 1$......0.2ccc0stececeecsee sve eect cere Volezbn Crile) 
lente ne x cistance to calculate Over 1S........,.......-<...---.-- %3.2f m\n',h) 
fprintf(‘The y distance to calculate over is................cccce. %3.2f m\n',h) 
imme ene number Of X prid POINts 1S.......:..0........0...0.. %4.0f\n',stop-1) 
fprintf(‘The number of y ae POGUES SS ieee te ss stse cs %4.0f\n',stop-1) 
fprintf(‘The x step size is.. eer. ole tn .dx) 
ION IMPNMnIGRy SLC) SIZE 1S..¢..00.2...6-2ebscistedes cess cee % 1.4f m\n',dx) 
fprintf('------------------------------------------------------------------ \n') 


Ip OR Sea cats cc eacr a acme nea AS RE oh a A HE Ae ae ee ke a he He Ne ek ae 


*\n') 


fprintf(’ Conditions ahead of the shock 

\n') 

Intel tnt Gc nT MIR lel Semen PaO IT... slelas ease. %1.1f\n',g) 
POMINL(euhe MttAnplesSULe 1S) 22 ..cs cesses. %3.3g Pa\n',Pin) 
Home iiemmltralGensity ($428.6 506<..4,.0000scstes %1.4f kg/m‘3\n' rho) 
TOME We MMETAlNPeIOeTAUUNG 1S tye cast... ..c.-1--. 2s. %3.3g K\n',T1) 
font Whe SpeediousOuUldis..a00.-............:.........763.41 m/s\n',c) 


Mapes ES TOE F e  I  A AE AK RK A ACC eK 2 Ok A 


*\n') 
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endt=input(‘Enter the number of time steps to calculate over: '); 


if endt<=0 
error(’You must enter a number of time steps larger than zero!') 
end 

ON ee 

% define vector size 

for k= -stop 

for y=1:stop 


op(x,y)=0.0; 
np(x,y)=0.0; 
Olr(x,y)=0.0; 
nr(x,y)=0.0; 
br(x,y)=0.0; 
Olru(x,y)=0.0; 
nru(x,y)=0.0; 
Olrv(x,y)=0.0; 
nrv(x,y)=0.0; 
bru(x,y)=0.0; 
nu(x,y)=0.0; 
brv(x,y)=0.0; 
nv(x,y)=0.0; 
oE(x,y)=0.0; 
nE(x,y)=0.0; 
bE(x,y)=0.0; 


end 
end 


% Calculate Peak static overpressure Ps based on charge 
% size and range to charge the driving force for the 
% calculations to follow 


% 1=range to charge 


1=input(Enter the range to the charge in meters: '); 
if r<=0 
error(You must enter a range larger than zero!’) 
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end 
% Print out table of explosive types 
% T=type of explosive 
fprintf('-------------------------------------------------------- \n') 
fprintf(’ Choose the explosive using the table below\t\n') 
Iefor AHA (GAG: @OTNPIO UGE 32 seslethsa sac sdinwsssacocdec cote <4t- l\n') 
fepemilin gan UDC ( COMICRORICC) oscars ctcacsvesetsca-sodssvicle cavers. 2\n') 
IyfevTnU ROT Wal UX ete ees wage tess: craessccdosetas.voceess--00es- 3\n') 
Pom NITmOP ty Ceti (GUI) i.i.c....0cs0cc+-es0s 0400-0202 A\n') 
Hyp TAUTMAAL (QMUMINIE rAd cos. Pox, Ge. ieala Nau acdsusidas ade-a4 0%: S\n') 
ionlumnemlastame Gelatin ....0.:510ceceees ieecedesseens 6\n') 
fprintf(‘60 percent Nitroglycerin dynamite ..................... 7\n') 
MP NPAT CS ROU KO ee oe es ane cc sacs deadesdacdicestiessass. 8\n') 
fprintf('-------------------------------------------------------- \n') 
Wor vrcpre net ee 
% Ask for what type of explosive to use 
i) rn epee nhs ere eee 
=input(‘enter the type of explosive to use: '); 
if T==1 
S=1.148; 
elseif 1—— 
S=1.185; 
elseif T== 
S=1.256; 
elseif T== 
S=1.481; 
elseif T==5 
S=1.000; 
elseif T= 
S=1.000; 
elseif T== 
S=0.600; 
elseif T== 
S=1.250; 
else 
error(’You must enter an explosive type!') 
end 


fprintf(‘The TNT conversion factor for this explosive is: %1.3f\n',S) 





ve 


% Ask for mass of charge 
(kg) 
0 a pe 
M=input(‘Enter the mass of the charge in kilograms: '); 
if M<O 
error('You must enter a mass larger than zero!') 
end 
W=MiS; 
i ero oe Seren oe oe nnn 
% Calculate peak static overpressure (B) 
Eee IN enantio ____-_-...--- 
z=1/(W“(1/3)); 
Ps=(6.7/(2%3)); 
if Ps<10 
Ps=((0.975/z)+(1.455/(z*2))+(5.85/(z*3)))-0.0 19; 
else Ps=(6.7/(z%3)); 
end 
if Ps<0.1 
Ps=Ps; 
end 
tu ater se eee 
% Convert Ps(the static overpressure) from Bar (Pa) 
% to Pascals from overpressure to absolute 
% pressureand apply conversion factor for 
% ground burst 
ire ere ee 
Ps=1.8*((Ps*100000)+Pin); 
De ee Ee 
% Print output of values on screen 
i eee pe ee eee eee Oe Se ----------- 
fprintf('------------------------------------------------------------------ \n') 
fprintf(’ Explosive values 
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fprintf(‘The mass of the charge in TNT equivalents is.............%4.3f kg\n',W) 
pint Whe staticroverpiessure of the charge is................... %1.3g Pa\n',Ps-Pin) 
fprintf('------------------------------------------------------------------ ') 

DF pet 

% Calculate intial velocity using Ps 
ee ei 


Mach=sart(((Ps/Pin)-1)*((g+1)/(2*g))+1); 
Vin=c*(2/(g+1)*(Mach-(1/Mach))); 


% Print output of values on screen 
fprintf(‘The initial wave velocity 18................:.0.....%4.4f m/s\n', Vin) 
fprintf(‘The shock wave Mach number ig....................0c008 %2.4f\n',Mach) 
fprintf('------------------------------------------------------------------ \n') 

I he ee cpceene eek en ----- 

% Calculate values behind the shock 
i NINE Bor ey tn ts 


tho2=rho/(1-(2/(g+1))*(1-(1/Mach’2))); 
T2=Ps/(rho2*(8314.51/29)); 
c2=(c*2*T2/T1)10.5; 


Ve eee cs ce act ee Ee A Re oe ok ee BE eK Seale 


in) 
fprintf(’ Conditions behind the shock 
\n') 
ioiahiatit MMe naTeSSUNe (Sey eee a.........Y0l. 3p Pa\n',Ps) 
MOMMCR( Mune GEMSILY 19...0.2.2....... MYe...s..cec-ceescdecasesses %1.4f kg/m‘’3\n',rho2) 
ipaeatat (CMM eG PEMMAPETAUINS (Siui,casiliee.sc.c.cvvsssesesvoutebaareves- %3.3g K\n',T2) 
MOU MMUUNE Mp SEC Ol SOUME 1S soys.c28 20sec cc ese¥ses.srsees ee %3.4f m/s\n',c2) 
HePOO ToL Tecra ck i RE EE AE Ree ok ok RR oR RE REE ERK EEE 
*\n') 
ea SRC 2 
i eee eee ee Mines eer ee 2 tne 
% Define initial conditions 
eee eeeemuen ewer Sey ___.-.-.--- 
for x=I:stop 
for y=l:stop 
op(x,y)=Pin,; 
Olr(x,y)=rho; 


Olru(x,y)=0; 





a 


Olrv(x,y)=0; 
oE(x,y)=Pin/g1; 
end 
end 
ee tee Oe a oo Se ee 
% Define boundary conditions at grid point x=1 y=1 


I 
Vee, 
x=2; 
op(x,y)=Ps; 
Olr(x,y)=rho/(1-(2/(g+1))*(1-(1/Mach’2))); 
Olru(x,y)=Olr(x,y)* Vin; 
Olrv(x,y)=Olr(x,y)* Vin; 
oE(x,y)=(Ps/g1)+(0.5*(Olru(x,y)*2+Olrv(x,y)*2)/... 


Olr(x,y)); 
DNR nee 2g ie oH bn gpa ar ecg acean acca 
% Calculate the new speed of sound based on new pressure 
% and density values 
DP ee 
nc=sqrt(g*op(x,y)/Olr(x,y)); 
ga 
% Calculate the initial time step 
yf ee 
dt=CFL*(1/(abs(Olru(1,1)/Olr(1, 1))/dx+abs(Olrv(1,1)/Olr(1,1))/... 
dy cien (bce Zemly dy 2) (0.5))); 
OTTER Mee Atal EITMe SUED 1S ay. .ccc.43c000hs)2-+000--20- %1.S5e s\n',dt) 
0 pa Orne 8 ck 
% Initialize total time counter 
0 an nn Nelda 
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stop=stop-1; 


for t=l:endt 


Vy---- nnn nn nnn nnn nn on nn nnn nnn = = 2-22-22 
% NBS for column 1| and row 1 
tauren re ee ee 
dtmin=1.; 
dtdx=dt/dx; 
dtdy=dt/dy; 
0 pe dig Peewee oe) ee 2k. 
% Calculate total time 
Wp---------------------------------~------------------------------- 
tt=tt+dt; 
Ofy---------------- 2-2-2 = -------------------------------------- 
for x=1:stop 
for y=1:stop 
eee ne ee eee ne. eT 
% PREDICTOR 
0 IE Ege een nnn 
% Solve as a function of time and position for 


% increasing time 
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% a ee 
of Density(r) 
% ee a Se tes ein as ets i ew ees 
br(x,y)=Olr(x,y)-... 
(dt*(((Olru(x+1,y)-Olru(x,y))/dx)+... 
((Olrv(x,y+1)-Olrv(x,y))/dy))); 
br(x,stop+ 1)=br(x,stop); 
br(stop+1,y)=br(stop,y); 
br(stopt+1,stopt+1)=br(stop,stop); 
0% ea the Se ae a re Bt Se 
0% Momentum(ru) 
% ee ee ee ae ee ee aaa oe eee be de daeeeakeiuem 
bru(x,y)=Olru(x,y)-dt*... 
(OH dx * 
(Olru(x+1,y)*2/Olr(x+ 1,y)-... 
Olru(x,y)*2 /Olr(x,y) +... 
(op(x+l,y)-op(x.y) —_)))+... 
(lidy == 
(Olrv(x,y+1)*Olru(x,y+1)/Olr(x,yt+1)-... 
Olrv(x,y) *Olru(x,y) /Olr(x,y)_))); 
bru(x,stopt+1)=bru(x,stop); 
bru(stopt1,y)=bru(stop,y); 
bru(stopt1,stop+1)=bru(stop,stop); 
% a ee Se ee eames 
of, Momentum(rv) 
% a Nc 


brv(x,y)=Olrv(x,y)-dt*... 
(lax 
(Olru(x+1,y)*Olrv(x+1,y)/Olr(x+1,y)-... 
Olru(x,y) *Olrv(x,y) /Olr(x,y) _))+... 
(1/dy*... 
(Olrv(x,y+1)*2/Oli(x,yt+1)-... 
Olrv(x,y)*2 /Olr(x,y)... 
(op(x,y+ I )-op(x,y) )))): 


brv(stop+1,y)=brv(stop,y); 
brv(x,stop+1)=brv(x,stop); 
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brv(stop+1 ,stop+1)=brv(stop,stop); 


OP pa 
% Energy(e) 
ee re ee ee eae Pee a Sea koe - 
bE(x,y)=o0E(x,y)-dt*... 
(Giiax*.; 
((oE(x+1,y)+op(xt+1,y))*Olru(x+1,y)/Olr(x+1,y) -... 
(oE(x,y)+op(x,y) )*Olru(x,y) /Olr(x,y) )+... 
1/dy*... 
((oE(x,y+1)+op(x,y+1))*Olrv(x,y+1)/Olr(x,yt+1) -... 
(oE(x,y)t+op(xy) )*Olrv(x,y) /Olr(x,y) ))); 
bE(stop+1,y)=bE(stop,y); 
bE(x,stopt 1)=bE(x,stop); 
bE(stop+1,stopt+1)=bE(stop,stop),; 
1 (Ns Ee ek _s 
% return for next x 
end 
% retum for next y 
end 
0 eae enna elt ODP a a oat en 
for x=2Z:stopt+| 
for Vee StOpial 
(7 REN a 
% CORRECTOR 
pe Ne 
of Density(r) 
nen NENe IMME 


nr(x,y)=(br(x,y)+Olr(x,y))/2-... 
( dtdx/2*(bru(x,y)-bru(x-1,y))+... 
dtdy/2*(brv(x,y)-brv(x,y-1))); 


nr(1,y)=nr(2,y); 
nr(x, 1)=nr(x,2); 
nr(2,2)=Olr(2,2); 
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nr(1,1)=Olr(2,2); 


% Left wall boundary conditions 


% Momentum(ru) 


Px=g1*(bE(x,y)-0.5*((bru(x,y)*2/br(x,y))+(brv(x,y)*2/br(x,y)))); 
Pxm=g1 *(bE(x-1,y)-0.5*((bru(x-1,y)*2/br(x-1,y))+(brv(x-1,y)*2/br(x-1,y)))); 


nru(x,y)=(bru(x,y)+Olru(x,y))/2-... 
((dtdx/2*... 
(( bru(x,y) *bru(x,y) /br(x,y) +Px )-... 
( bru(x-1,y)*bru(x-1,y)/br(x-1,y)+Pxm)))+... 
(dtdy/2*... 
( bru(x,y) *brv(x,y) /br(x,y)-... 
bru(x,y-l)*brv(x,y-l)/br(x,y-1) _))); 


ee a ae a ee 
nru(1,y)=nru(2,y); 
nru(x, | )=nru(x,2); 
nru(2,2)=Olru(2,2); 
nru(1,1)=Olru(2,2); 
% Momentum(rv) 
0 ie ee ee ee ee eee eee ee = o-oo n= 


Py=gl*(bE(x,y)-0.5*((bru(x,y)*2/br(x,y))+(brv(x.y)*2/br(x,y)))) 
Pym=g 1 *(bE(x,y-1)-0.5*((bru(x,y-1)*2/br(x,y-1))+(brv(x,y-1)*2/br(x,y-1)))); 


nrv(x,y)=(brv(x,y)+Olrv(x,y))/2-... 
((dtdx/2*... 
( bru(x,y) *brv(x,y) /br(x,y) — -... 
bru(x-1,y)*brv(x- 1,y)/br(x-1,y) ))+... 
(didy/2™... 
( (brv(x,y)*brv(x,y)/br(x,y)+ = Py )-... 
(brv(x,y-1)*brv(x,y-1)/br(x,y-1)+Pym)))); 
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nrv(1,y)=nrv(2,y); 
nrv(x, |)=nrv(x,2); 
nrv(2,2)=Olrv(2,2); 
nrv(1,1)=Olrv(2,2); 


% Left wall boundary conditions 


0% Energy(E) 


nE(x,y)=(bE(x,y)+oE(x,y))/2-... 
(iia 273: 
((bru(x,y) /br(x,y) *(bE(x,y) +Px )) -... 
(bru(x-1,y)/br(x-1,y)*(bE(x-1,y)+Pxm)))) +... 
(didy2= 
((brv(x,y) /br(x,y) *(DE(xy)+Py )) -... 
(brv(x,y-1)/br(x,y-1)*(bE(x,y-1)+Pym)))) _ ); 


Yl 208 Pele ere ee ek 
se cen aoe SN 
nE(1 sy)=nE(2,y); 
nE(x,1)=nE(x,2); 
nE(2,2)=oE(2,2); 
nE(1,1)=o0E(2,2); 
OAs eon era so ope een teneren ES NN 
“return for next x 
end 
07ers meneame entemneek 
%return for next y 
end 
ene 0 er 2 ne 
for y=I:stopt+l 
for Sey ae 
jn 8 eee 
% Calculate pressure using the continuity equation for 
% the ideal gas law and the new values for energy, 
% momentum and density 
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np(x,y)=(g1)*(nE(x,y)-(.5*((nru(xy)2/... 
nr(x,y))+(nrv(x,y)*2/nr(x,y))))); 


nN ete ie ew atamen heels 20000 0g 8 7 ioe 
% Calculate velocity using the new values for momentum 
% and density 

(EI aot SAE AD cae a 0a sos 


nu(x,y)=nru(x,y)/nr(x,y); 
nv(x,y)=nrv(x,y)/nr(x,y); 


Ue a aa a a a es asd ened tan 


% Check stability 
% Calculate the new speed of sound 


dtn=CFL*(1/(abs(nu(x,y))/dx+abs(nv(x,y))/dyt... 
nce*(1/dx*2+1/dy*2)*(0.5))); 
if dtn < dtmin 


dtmin=dtn; 

% end dtmin if statement 

end 
Tinie f nc eeee ees ee ee 
a a nT OS se 
5 aaa nN Lo TE 
% return for next x 

end 
5 an RCE en nce 


“return for next y 
end 
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% Update time step and CFL(Courant number) 


ja eaenaeao—saao==- Pp pg Ea SE pe IE EE i a 
dt=dtmin; 


fprintf(‘delta t= %1.S5e s\tTotal time= %1.5e s\n’,dt,tt) 


(eed. oc cee ses oe seed oe oe on. ---2- 
% Output the values for pressure, density and velocity 
% using graphs 
fe 
% Pressure 
for Ke StOpa! 
for y=2:stoptl 
press(x-1,y-1)=np(x,y); 
end 
end 
colormap(cool) 
surf(press) 
hold on 
pcolor((press)-Pin) 
*%contour(press,20) 
title(’‘Absolute Pressure’) 
xlabel('X grid points’) 
ylabel(‘Y grid points’) 
Zlabel('Pascals') 
pause(.01) 
%hold off 
De 
% Replace old values with new values and begin 
% next time step 
0 un ee 2 5 enn ee 
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Olr=nr; 
Olru=nru; 
Olrv=nrv; 
oE=nE; 
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